A Symplectic Integrator for non-integrable Hamiltonian Relativistic Systems 
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1. INTRODUCTION 

In a few years from now, a new era of astronomy is 
expected to begin when the Advanced LIGO and other 
gravitational wave detectors are anticipated to receive 
the first signals. The main source of gravitational wave 
signals are binary systems consisting of compact objects 
whose mass proportion ranges from equal masses to ex- 
treme mass ratio. The gravitational wave signals from 
binary systems are closely related to the motion of the in- 
spiraling compact objects. Thus, orbital studies of binary 
systems have attracted much attention (e.g., [1, 2]). In 
general, orbital motion can be studied by using a Hamil- 
tonian function which splits into an orbital and a spin 
part. 

One possible approximation is the post-Newtonian 
(e.g., [3] and references therein). In the post-Newtonian 
approximation the Hamiltonian's orbital part is ex- 
panded in powers of -7. This approach has also been 
extended to binaries that are perturbed by a third, 
much lighter object (e.g., [4]). These post-Newtonian 
approaches have in common that the mass ratio of the 
binary should not be very large. 

When the binary system consists of compact objects of 
extreme mass ratio, the motion in such a system is called 
an Extreme Mass Ratio Inspiral (EMRI). We expect to 
find such binaries at the center of galaxies, where a com- 
parably light compact object inspirals into a central Su- 
permassive Black Hole (SMBH). The lighter inspiraling 
compact object moves adiabatically from one geodesic 
orbit to another. Thus, an EMRI can be approximated 
by the geodesic motion of a test particle in the spacetime 
background of the central supermassive object (e.g., [5] 
and references therein). 

One interesting aspect of EMRI is that the shape of 
the spacetime background is encoded in the gravitational 



waves emitted by the inspiraling object. Ryan showed 
that we can, in principle, extract this information from 
the gravitational waves [6]. Therefore, when future low 
frequency space-borne gravitational wave detectors are 
launched, we will be able to test whether the spacetime 
backgrounds around SMBHs are indeed described by the 
Kerr metric or by something else [5] . In order to provide 
such a test, Collins and Hughes [7] constructed a per- 
turbed Schwarzschild spacetime background and called it 
a bumpy black hole. Since then, several authors have sug- 
gested different tests [8-21] using various forms of bumpy 
black holes (e.g., [22-24]), which are also known as non- 
Kerr metrics. 

Usually, these non-Kerr metrics are stationary and 
axisymmetric, but they miss the fourth integral which 
would make them integrable systems, as the Carter con- 
stant [25] does in the case of the Kerr metric. Thus, 
in general the non-Kerr metrics correspond to non- 
integrable systems, which in practical calculations of- 
ten are investigated for chaos (e.g., [14, 26-29]). This, 
in turn, implies numerical simulations over long time- 
intervals. Therefore, high-efficient integrators with a 
good long-time behavior are required. 

Over the last decades, part of the Numerical Analy- 
sis community has focused on long-term integrations of 
differential equations. Geometric or structure preserv- 
ing algorithms, such as symplectic methods for Hamilto- 
nian equations of motions or symmetric integrators for 
reversible systems, have been developed aplenty. Con- 
trary to standard explicit integrators as Runge-Kutta 
schemes, these integrators nearly conserve first integrals 
(e.g., the energy) and overall integration errors increase 
only slowly with time. Whereas for standard integration 
schemes, the overall error is normally proportional to the 
square of the length of the integration interval t% , it only 
increases linearly with ti for structure preserving integra- 
tors. A detailed discussion of such methods can be found 
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The efficiency of numerical integrators can be increased 
by the use of adaptive step size controllers. But, by 
using standard step size controllers for geometric algo- 
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rithms, we lose the integrators' good long-time behav- 
ior. However, over the last 15 years, algorithms allowing 
structure preserving step size control have been devel- 
oped ([31, 32]). This development has resulted in a lot of 
available tools to confront problems in classical astron- 
omy. 

However, in general relativity, the structure of the 
equations is more complex. The inclusion of spins to 
a binary system leads to a Poisson-system. For the 
post-Newtonian formulation of such a system, a non- 
canonically symplectic integrator has been found, sec 
[33]. In the geodesic approach, the Hamiltonian of the 
orbit has a position dependent mass-matrix. Thus, from 
this point of view, the system of differential equations 
is non-separable. Furthermore, the vector field of the 
resulting equations of motion is strongly varying in non- 
integrable spacetime backgrounds. These two features 
pose a new set of problems for the implementation of ef- 
ficient and accurate algorithms. As an answer to these 
problems we introduce a new integration scheme for such 
differential equations based on Gauss-Runge-Kutta col- 
location methods. In particular, we use a new struc- 
ture preserving adaptive step size controller. As a re- 
sult, we get a symmetric, reversible integrator which ef- 
ficiently conserves the constants of motion. We tested 
this integrator's performance in the case of the Manko, 
Sanabria- Gomez, Manko (MSM) metric [24] and the re- 
sults show that the new integration scheme is much faster 
than an adaptive step fifth-order Cash-Karp Rungc- 
Kutta scheme. 

The paper is organized as follows: Section 2 summa- 
rizes some basic elements about the MSM metric, the 
geodesic motion in the MSM spacetime background. In 
the subsequent Section 3, we consider the topic from a 
Numerical Analysis point of view and discuss the char- 
acteristics of the geodesic equations. Our step size con- 
troller, along with the integrator, is presented in section 
4. The performance of the integrator in numerical tests 
is shown in section 5. Section 6 summarizes the main 
results of the study. 



2. MANKO, SANABRI-GOMEZ, MANKO 
SPACETIME BACKGROUND 

2.1. The metric 

We tested the various integration schemes on the five- 
parameter vacuum solution introduced by Manko et al. 
[24]. The MSM solution is asymptotically flat, axisym- 
metric and stationary, and describes the "exterior field 
of a charged, magnetized, spinning deformed mass" [24]. 
The MSM spacetime depends on five real parameters, 
namely the mass to, the spin (per unit mass) a, the total 
charge q, the magnetic dipole moment Ai 1 and the mass- 
quadrupole moment Q. M. and Q are represented in the 



metric by the real parameters p, and 6, i.e. 

M = p + q(a — b) , 
Q = -m(d-S -ab + a 2 ) , (1) 

where 

p 2 — m 2 b 2 
m 2 — (a — b) 2 — q 2 

d := \[m 2 -{a-b) 2 -q 2 ] . (2) 

The metric for the MSM spacetime can be given by 
the Weyl-Papapetrou line element 

ds 2 = -f(dt-cud(f>) 2 

+ r 1 [e 2 i(dp 2 + dz 2 )+p 2 dc}> 2 ] , (3) 

where all metric functions /, uj, and 7 are considered 
as functions of the prolate spheroidal coordinates u,v, 
while the coordinates p,z are the corresponding cylin- 
drical coordinates. The transformation between the two 
coordinate systems is 

p = K\J (u 2 — 1)(1 — v 2 ), z = kuv , (4) 

where 

k := Vd+~5 . (5) 

The metric functions are 

/ = E/D , 

e 2 '< = E/16k 8 (u 2 - v 2 f , (6) 

uj = (v 2 - l)F/E , 

where 

E = R 2 + A1A2 5* 2 , 

D = E + RP+X 2 ST , (7) 

F = RT-XiSP , 

Ai=k 2 (u 2 -1), \ 2 =v 2 -l (8) 

and 

P := 2{Kmu[(2KU + m) 2 -2v 2 {25 + ab-b 2 ) 

- a 2 + b 2 -q 2 ]-2 K 2 q 2 u 2 -2v 2 (A5d-m 2 b 2 )} , 
R := 4[k 2 (u 2 -l) + S(l-v 2 )] 2 

+ (a - b)[(a - b)(d - S) - m 2 b + q p](l - v 2 ) 2 , 

(9) 

S := -4(a-b)[ K 2 (u 2 -v 2 ) + 2Sv 2 ]+v 2 (m 2 b-q p) , 
T := A(2Kmbu + 2m 2 b-q p)[n 2 (u 2 -l) + 5{l-v 2 )] 
+ (1 -v 2 ){{a-b){m 2 b 2 -ASd) 

- {Anmu + 2to 2 - q 2 )[{a - b)(d-5) - m 2 b + q p]}. 

The MSM metric is of astrophysical interest as it has 
been suggested to be appropriate to model neutron stars 
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[24, 37] . By investigating the dynamics of the geodesic or- 
bits in the MSM spacetime, Dubcibe et al. [26] and Han 
[27] debated on the appearance and on the significance 
of chaos in the MSM model. Even if the main purpose 
of our article is to present a new integration scheme, in 
section 5 we contribute to the aforementioned debate. 



As the system is autonomous (^r = = 0), the Hamil- 
tonian function is an integral of motion and equal to 
H = -mo/2. 

3. STANDARD SYMPLECTIC SCHEMES 



2.2. Geodesic motion 

The equations of geodesic motion of a "test" particle 
of rest mass Too in a spacetime given by the metric g^ u 
are produced by the Lagrangian function 

L = i to g^u x>*x v , (10) 

where the dot denotes derivation with respect to proper 
time t. The Lagrangian (10) has a constant value L = 
—m /2 along a geodesic orbit, due to the four-velocity 
9tiv x^x v = — 1 constraint. 

Since the spacetime is axisymmetric and stationary, 
the corresponding momenta 

are conserved. These are the specific energy 

1 dL 
m di 

and the specific azimuthal component of the angular mo- 
mentum 



E 



(12) 



L - — — 

z toq dip 



(13) 



For brevity, we refer hereafter to these two integrals sim- 
ply as the energy E and the angular momentum L z . Due 
to these two integrals of motion, we can restrict our study 
to the meridian plane. We just have to re-express (12), 
(13) in order to get i and <j> as functions of E and L z , 
and then replace i, <fi in the two remaining equations of 
motion. Then, from the original set of 4 coupled second 
order ordinary differential equations (ODEs), we arrive 
to a set of 2 coupled ODEs. 

By simply applying the Legendre transform 



H = p M i" - L 



(14) 



to the Lagrangian function (10), we get the Hamiltonian 
function 



H = 



2toq 



(15) 



where the momenta p^ are given by (11) and p v — mox v . 
The Hamiltonian equations are 

,.dH 



dp" 



(16) 



,dH 
dx" 



We have seen in section 2.2 that the geodesic equations 
of motion can be described by a Hamiltonian formalism. 
Thus, symplectic schemes should be appropriate for in- 
tegrating these equations. 



3.1. Notation 

In the rest of the paper we drop the Tensor Analysis 
covariant and contravariant symbolism by indices and ex- 
ponents, and we adopt a symbolism more convenient for 
Numerical Analysis. By bold characters, we denote vec- 
tors, and, by arrows over characters, we denote vectors 
of vectors. 

We define 



y := 



I 



1 " 1 -I II. 

/(y) := J-'VHiy) 
and write the given Hamiltonian system as 

%-tW. (17) 

Further, $^ denotes an integration step of step size 
h. If the current position in phase space is y„, then 
propagates the system to the next position y„+i, i.e. 

y™+i = $/.(y n ) ■ 

To simplify the Runge-Kutta-schemes notation we define 
the vectors Z i; the functions /(Zj) and the coefficient 
matrix A by 

Z:=(Zj ... Z S ) T , 
F(Z):=(f(y n + Z 1 ) ... f(y n + Z s )) T , (18) 

(an ... a u 
: : 
a s \ ... a s . 

The exclamation mark "!" over relation symbols de- 
notes requirement. Furthermore, X denotes the phase 
space and X denotes the space on which the Z are de- 
fined, i.e. X := X x ... x X. 

As a norm for matrices we use the Frobenius norm, 
which for a matrix A is 



* ,3 
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3.2. The equations of motion from a numerical 
point of view 

From a numerical point of view, it is easier to handle 
the Hamiltonian equations (16) in the form (17). The 
Hamiltonian system (15) is non-separable in the sense 
that the mass-matrix in the Hamiltonian depends on the 
positions x and, therefore, the momenta cannot be sep- 
arated from the position coordinates. In the case of a 
non-integrable system, the set of coupled first order dif- 
ferential equations (17) can be very sensitive to changes 
in the arguments y. This happens when we evolve chaotic 
orbits by In other words, if Df(y) denotes the Ja- 
cobian of /(y) at a point y of the phase-space X, then 
there exists a subset U C X, so that 

||£>/(y)||»l VyeW . (19) 

As the Hamiltonian system (15) is not separable, a 
natural choice for a symmetric, symplectic integrator is a 
Gauss- collocation scheme. An s-stage Gauss-collocation- 
scheme is an s-stage Runge-Kutta scheme 

s 

y n +i = y n + h^2bif(y n + Z;) , 

1=1 

s 

Z i = h^a i jf{y n + Z j ) , (20) 
whose inner stages c, are given as 

Ci = Tjl 1 + c ») . 

where q are the roots of the Legendre-polynomial of de- 
gree s. With the help of the Lagrange-polynomials 

Cj Cj 

the weights of the scheme are given as 




Using the notation introduced in eqs. (18), the system 
of implicit eqs. (20) is rewritten in a shorter form 

Z = h{A®I)F{Z) , (21) 

where <g> denotes the direct product of matrices. 

There are two possibilities to solve this implicit equa- 
tion for the inner stage values Z. 

• The simplest and the most popular method is the 
Fixed-Point iteration. One has to solve iteratively 

Z k+1 = h(A®I)F(Z k ) . (22) 



The Banach fixed-point theorem guarantees that 
the iteration convergences towards the correct Z, 
if it is a contraction. This requires 

\\z k+2 - z k+1 \\ < \\z k+1 - z k \\ . 

Because of 

\\Z k+2 - Z k+1 \\ = \\h(A® I){F(Z k+1 ) - F(Z k ))\\ 
<max \\h(A<E)I)DF(Z)\\\\Z k+1 - Z k \\ , 

this implies 

h\\(A® I)DF{Z)\\ < 1, yZeX. (23) 

• If the convergence of the fixed-point iteration is not 
ensured, we can look forward to get better accuracy 
by searching for the roots of 

F(Z) := Z -h(A®I)F(Z) (24) 

via a modified Newton iteration. In this case, we 
have to iterate 

Z k+1 = Z k + AZ k (25) 

AZ k = -M- 1 F(Z k ) , (26) 

with 

M:=I-(A®I)(I®Df(y n )) . (27) 

Again, convergence can only be ensured for a con- 
traction. Therefore, we need 

\\Z k+2 - Z k+1 \\ = 

\\Z k+1 - M- 1 F(Z k+1 ) - (2 k - M- x F{Z k )} || 

<max h\\I - M- x DF{Z)\\\\Z k + x - Z k \\ 
zex 

< \\Z k+1 -Z k \\ , 
or 

hWl-M^DF^W < 1, VZ e X . (28) 

As a single iteration step for the fixed-point-iteration 
is computationally much cheaper than in the newton- 
iteration case, one should prefer the first whenever it is 
applicable. 

In both cases, inequality (19) together with the respec- 
tive requirements (23) or (28) lead to severe restrictions 
on the possible step-size h. 
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3.3. Existing step-size controls 

For efficiency's sake one would expect to use small time 
steps only for the subspace for which condition (19) holds. 
Therefore, a reliable step size control is necessary. A 
feasible step size control has to preserve the geometric 
properties of the underlying integration scheme. As it 
was shown in [34] , there cannot be an efficient symplectic 
integrator with adaptive step size. Thus, other structure 
preserving properties of the integrator like symmetry, i.e. 



<S>_ h (y)=<S>-\y) 



and reversibility, i.e. 



<P h 1 o C = C ° ®h 



(29) 



(30) 



for the involution C(P> X ) = ( — P> x )> have to be con- 
served. If we apply such an integrator to a symmetric, 
reversible differential equation, we would expect the same 
advantageous long-time behavior as for a symplectic in- 
tegrator, e.g., [30], chapter VIII. 

One of the most efficient step size controls has been 
proposed by Hairer and Soderlind [32]. They express 
the variable step size sequence r = n, r„, r n +i, ... 
by means of a constant step size sequence e = 
6\, e n , e n +i, ... and a scalar function er(y) via dr — 
a(y)de. For the dependence of y on e, they obtain 



dy _ dydr _ dr 
de ~ dr de ~ J{y 'de 



Defining 



one finds 



w := 



*(y) 



dw 1 dy 

de a 2 de 

= --V<7/(y) 



(31) 
(32) 

(33) 
(34) 



(35) 



We now combine eqs. (31) and (33) and get the system 

{-^af(y) =: G(y) / 

Solving this system with the scheme 

w n+ i = w n + e-G(y n ) (36) 
y„ + 1 = $_^(y„) (37) 

w n + i =w n+ i +e-G(y„+i) , (38) 
1 



Wq = 



o-(yo) 



(39) 



we get a symmetric, reversible integrator, as it is shown 
in [32]. 



The above integration scheme is explicit in the variable 
w. Therefore, eq. (36) represents a constraint on possible 
step size adapters c(y). If we don't want the underlying 
step size e to be too small, G(a(y)) has to be bounded. 
This is similar to the case when we use explicit Runge- 
Kutta schemes for stiff differential equations, where, due 
to a very small range of stability, the step size has to be 
unfeasibly small. 

However, for the iterations to converge, we have to 
choose a step size small enough to get a contraction. 

• If we want to solve equation (21) with a fixed-point- 
Iteration, the step size must be somehow propor- 
tional to the inverse of ||D/(y)||, i.e. we have 



1 



a(y) cx 



IA/(y)ll 



(40) 



and such a controller a yields 



, , E«> (Df(y)) v (^i^ 1 ) (/(y»* , > 

GW ■ 

For regions of the phase space, where 

||£> 2 /-/ll»llAfll , (42) 

this would conflict with the requirement which de- 
mands G(y) not to be too large. Numerical tests 
even show that in such cases the step size can be- 
come negative. 

• In order to guarantee the convergence of the 
newton-iteration, i.e. to satisfy condition (28), it 
is required that 



cr(y) cx 



1 



[7- M- 1 DF(Z)|| 



||7 - (7 - (A ® D/(y„)))"i(7 - (A ® I)DF{Z))\\ 

(43) 

It is by no means clear how one can consider the 
last expression as a function of y. The only pos- 
sibility would be to set Z = (0 . . . 0) T in the last 
expression. The calculation of the corresponding 
G(y) would then lead to an expression with a fac- 
tor 

1 

\I - (I - (A i £>/(y n )))-i(/ - (A i I)DF{Z)W ' 

Because of DF(Z) = I ® Df(y) for the Z chosen 
above, this factor would be 

1 



[7 - (7 - (A ® D/(y„)))-i(7 - (A ® Df(y)))\\< 



(44) 
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Thus, the evaluation of the function G at the point 
y„ (which is necessary in eq. (36)) is impossible as 
we would have to divide by zero due to the factor 
(44). We see that it is not possible to use the modi- 
fied newton iteration along with the presented step 
size control algorithm. 



because the collocation polynomial for y n +i = &h(yn), 
which is the interpolation polynomial through the points 
(ci,y ra + Zj), is the same as the collocation polynomial 
for y n = $_/j(y n+ i), i.e. the interpolation polynomial 
through the points (cj, y n +i + Zj). Due to property (49), 
we get 



4. A NEW ADAPTIVE STEP SIZE CONTROL 
OF MOTION 

We now present a new integration scheme $^ with an 
adaptive step size h(e,y). The integrator uses an s-stage 
Gauss-collocation method as the underlying integrator. 
The step size depends on an underlying constant step size 
e and the actual state of the system y. The integration 
scheme has the following properties: 

• For the step size h an equation similar to 

1 



h(e,y) oc 



\\Df(y)\\ 

holds, and, thus, condition (23) is satisfied. 



(45) 



• The integrator is both symmetric and reversible. 
Hence, the theoretical results on the long-time be- 
havior of such integrators are applicable (e.g., [30], 
chapter XI). 

The main idea of the step size control algorithm is to 
replace requirement (23) by 



|fc--D/(y)ll 



(46) 



for some e < 1. We slightly modify this formulation to 
\\^(Df(y n + Z 1 ) + Df(y n + Z s ))\\=e , (47) 

with Zi( s ) being the inner-stage value at stage l(s) of 
the Gauss-Runge-Kutta scheme. In the subsequent para- 
graphs we demonstrate that this modification gives the 
required structure preserving properties. 

For the integrator to be symmetric, we must make sure 
that the step size when propagating the system back- 
wards in time (i.e. from y n +i to y„) has the same be- 
havior as when propagating forward (i.e. from y„ to 
y„+i), i.e. we need 



h(e,y„) = -h(-e,y n+ i) 



(48) 



This, together with the symmetry of the underlying 
integrator (i.e. the validity of property (29) for constant 
steps h), ensures the symmetry of the whole integration 
scheme. Let Zj denote the inner-stage values for the inte- 
gration backwards in time. The symmetry of the Gauss- 
collocation method then results in 



Zj + y„+i = y„ + Z s+ i_j Mi = 1, 



(49) 



h(-e,y n+ i) 



!i[i?/(y„ +1 + Z 1 ) + A/(y« + i + Z s )]|| 

— e 

= mDf(y n + Z s ) + Df(y n + Z 1 )}\\ 

= -h(e,y n ) . (50) 

If we consider the new variable step size integrator 
<J>/j(e,y) as a constant step size integration scheme 
then the reversibility condition (30) reads 



*7 1 oC = C°*e • 



If we have 



h(-e,( y„+i) = -h{e,y n ) 



(51) 



(52) 



then the reversibility of the method is a direct conse- 
quence of the reversibility of the underlying integrator 
(30). In order to prove condition (52), we denote by Z 
the inner stage values for the integration that starts at 
C Yn+i- We then notice, that the reversibility of the 
Gauss-Runge-Kutta scheme reads 



Cy„+i + Zj = ((y n + Z s+ i_ 4 ) . 

This yields 

— a 

M-e,Cyn+i) = jj 



(53) 



{Df((y n+1 + Zi) + £>/(Cy„+i + Z s )]|| 
— e 

= ||ic[^/(yn + z s ) + A/(y™ + z 1 )]|| 

= -h(e,y n ) , (54) 

where the last equality is a consequence of the orthogo- 
nality of the involution £. 

Equation (47), together with equation (21) for the 
inner-stage values Zj, leads to the system of equations 



h(A®I)F(Z) 

mDf(Z 1 )+Df(Z s )]\\ 



(55) 



Inserting the second part into the first, we can easily 
apply a fixcd-point-Iteration, which yields 



Z k+1 



e{A®I)F{Z k ) 
£||D/(Zf)+£>/(Z*) 



(56) 
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FIG. 1: The relative error AH as function of the proper time r in semi-logarithmic scale. 



5. NUMERICAL EXPERIMENTS 

In this section we test our Integrator for Geodesic 
Equations of Motion (IGEM) and compare it with several 
other integration schemes. These are 

• A standard fifth order Cash-Karp-Runge-Kutta 
scheme with constant step size e (RK5con) as pro- 
posed by Numerical Recipes [35]. 

• The Cash-Karp-Runge-Kutta scheme (RK5var) 
with initial step size e and a step size algorithm 
as follows: 

- Starting with y„ calculate y„ + i = $h,(y„) 

— Calculate the Hamiltonian's relative error 



AH = 



H(y n+1 )-H 



H 



* If AH is greater than a minimum toler- 
ance value toll = 10 -12 , then repeat the 
integration step with h new — ^. 

* If AH is smaller than a second tolerance 
value to?2 = 10~ 14 , then double the step 
size h before calculating the next step. 

* If toli < AH < toh, then calculate the 
next step. 



• A standard symplectic integration scheme used in 
classical celestial mechanics given by the Gauss- 
collocation as underlying integrator together with 
the step size control explained in section 3.3 with 

°W = TOT ( CCM )' 



• Again the standard scheme from classical celestial 
mechanics with <r(y) replaced by cr(y) — 
(CCM2). 



IIAf(y)|| 



In all implicit schemes, we use a fixed-point iteration 
which we implement as proposed by [36]: In order to 
increase the accurary, one iterates until either the dif- 
ference between two subsequent Z is numerically or 
fluctuations due to round-off errors start to occur. In 
other words, we use the stopping criterion 



\Z 



fc+i 



or \\Z 



fc+i 



z*\\ > \\z k 



z 



fc-ll 



As test case for the algorithms, we take the system 
described in section 2. In particular, we use the set of 
MSM parameters and initial condition discussed in [27]. 
Accordingly, we choose a central object whose mass is 
m = 2.904, the spin is a = 1.549, the charge is q = 0, 
while the real parameters \x and b correlated with the 
magnetic dipole and the mass-quadrapole deformation 
are set to be and 0.8 respectively. Regarding the initial 
conditions, the energy and the angular momentum are 
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set to be E = 0.971 and L z = 9.3, respectively, while we 
fix p = z = 0.0. The integrators are tested for 3 different 
values of p, which correspond to 3 different cases of orbits. 
We propagate the system in time until the proper time 
reaches r/ = 500000 and compare the cpu times T ca i c the 
calculation lasted. As a measure of our code's accuracy, 
we track the relative error AH. We abort the simulations 
if AH > 10~ 6 . 



5.1. Case with p = 30.7 
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The initial condition p = 30.7 corresponds to a regular 
orbit. Fig. 1 shows that the error is much larger for the 
explicit schemes, although we have chosen a much smaller 
step size e. In table I we see that not only are the explicit 
integrators less accurate, but they also require much more 
computation time. Moreover, if we increase the step sizes 
e of the explicit schemes, the energy drift becomes much 
worse. On the other hand, the CCM is the fastest al- 
gorithm in this case because the corresponding regular 
orbit does not posses the special feature discussed in sec- 
tion 3.3 and expressed by the condition (19). Overall, the 
explicit integrators show a linear drift whereas the sym- 
metric integrators preserve the Hamiltonian's constant 
value. 



Integrator 
RK5con 
RK5var 
CCM 
CCM2 
IGEM 



c 

0.01 
0.01 
1.0 
1.0 
1.0 



Tcalc M 

222.9 
160.6 
17.3 
46.1 
41.3 



FIG. 2: The step size h of the IGEM as function of proper 
time r 6 [0,500]. 



the other hand, IGEM shows the best long-time behavior 
for affordable computational cost. 



CCM2 
CCM - 
RK5con 
RK5var 
IGEM 




100000 200000 300000 400000 



TABLE I: The cpu calculation times for the proper time in- 
terval r £ [0, 500000] for different integration schemes. 



FIG. 3: The relative error AH as function of the proper time 
r in normal scale. 



5.2. Case with p = 1.7 

Even though this initial condition (p — 1.7) corre- 
sponds to a regular orbit, the system of differential equa- 
tions behaves in a more ill-mannered way than in the 
p = 30.7 regular case. This is tshown by tracking the 
step sizes of the IGEM-integrator along the propagation 
(Fig. 2), as we have h cx i\ Df (l( T ) } \ \ - ||D/(y(r))|| is vary- 
ing fast with time. 

When we tested the RK5var scheme, it needed 1 week 
of calculation time only to arrive at r = 758. Therefore, 
we relaxed the restrictions on the acceptable relative er- 
rors per step and set tol\ — 10~ 9 and toh = 10 -11 . 
Now, if we compare the five integrations schemes in ta- 
ble II and Fig. 3 we notice that the symplectic standard 
schemes from classical celestial mechanics -although fast- 
show very bad conservation properties for the constant 
of motion. The relative error is even worse than for the 
constant step size RK5con which shows a linear drift. On 



Integrator e T ca i c [s] 

RK5con 10 -4 23813.2 

RK5var 10 -4 19207.7 

CCM 0.1 266.9 

COM2 0.1 2100.9 

IGEM 0.1 2202.1 

TABLE II: The calculation times Tcaic for the proper time 
interval r £ [0, 500000] for the different integrators. 



5.3. Case with p = 0.7 

The initial conditions with p — 0.7 correspond to a 
strongly chaotic region. The calculations become much 
more computationally expensive than in the other two 
previous cases. Therefore, we restrict the propagation 
interval to r £ [0,50000]. The runs using the well es- 
tablished integrators were aborted very soon due to low 
accuracy, or in the case of CCM because the variable step 
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size indeed became negative (table 77/)! 

The only integrator to pass this test was IGEM. For 
IGEM, wc plot the relative error AH in Fig. 4 and ob- 
serve that some peaks appear. This happens because, 
during each step, we have to solve the iteration (56). For 
this to be a contraction, we need 



D 



e(A(g>I)F(Z) 
UDf{Zx) + Df(Z a 



< 1 



(57) 



with 



D (•••); 



e{A®I)DF(Z)) 

/ i 



|p/(zi: 



(e(A« 



+ Df(Z 

I)F(Z) 



(58) 



i |A/(Zi) + A/(z s )ll 3 

x^(7J>/(Z 1 ) H +£»/(Z 5 ) H ) 



ap/(Zi) fci + A/(z s ) fc Q 

6% 



(59) 



The norm of the first term is smaller than 1 by con- 
struction of the algorithm but the second term can be 
large. This can cause non-convergence of the fixed-point 
iteration. As a consequence, the calculated orbit jumps 
to another trajectory. This can be visualized by plotting 
V 2 (t) := z 2 + p 2 in Fig. 5. We observe some peaks in 
the velocity. The calculation seems to loose its trajec- 
tory, erring around. But, due to the structure preserving 
feature of the integrator, the error in the Hamiltonian re- 
covers its previous value. Thus, the calculated orbit gets 
back to a nearby trajectory where it remains until erring 
for another short time. In total, only during a time range 
of 3 • 10~ 3 , out of the whole interval of 50000, the relative 
error is above AH — 10~ 9 and, thus, the algorithm con- 
serves the Hamiltonian for almost every r in the proper 
time interval. 

We tried to erase the peaks by using a newton-iteration 
instead of a fixed-point iteration. However, we observed 
that some peaks remained. Because the newton-iteration 
in this case requires the Hessian of /(y), the calculation 
time increased considerably to T ca i c = 4751.5 s. 

Our last numerical analysis argument is that IGEM is 
more convenient in computing Poincare sections than ex- 
plicit Runge-Kutta schemes. Without loss of generality, 
we suppose that, during the propagation y n — > y„+i, wc 
have crossed the z — plane. If we use an integration 
scheme of order p with step size s, we have y n and y n +i 
0(h p ) close to the real value. On the other hand, we get 
by Taylor-expansion 

y n +i = y n + 0(h) . 

Thus, whatever the order of the integration scheme, we 
only can give the location of the section with an accu- 
racy of 0(h), because we only know that it is somewhere 



3f 1e-0 



10000 20000 



40000 50000 



FIG. 4: The relative error AH for the IGEM method as a 
function of proper time r, semi- logarithmic scale. 




FIG. 5: The Velocity' V(r) 2 as function of proper time r. 



between y„ and y n +i. If we want to locate the sec- 
tion more accurately we can only go back to the point 
y„ and start the integration again with a much smaller 
step size. This, in turn, yields much more computational 
cost. Using the IGEM instead, we can take advantage 
of the collocation-property: It is known from Numerical 
Analysis that the solution at point y n +i calculated with 
an s-stage Gauss-collocation method coincides with 1(h), 
where 1(t) is the interpolation polynomial through the 
points (0,y„) and (c\,y n + Zi)...(c s ,y„ + Z s ). This is 
a special feature of collocation methods, which docs not 
hold for other Runge-Kutta schemes. Numerical Analy- 
sis shows further that the interpolation-polynomial 1(t) 
stays 0(h s ) close to the real solution within the whole 
interval between y„ and y„+i. Thus, if we search for 
the root of 1(t) (or rather the root of the z-component of 
1(t)) with a fast bisection method, we get the location of 
the Poincare section 0(h s ) close to the real section. Wc 
use this advantage of the IGEM to enter the discussion 
about the appearance of chaos in the MSM model. 

In [27], Han showed that chaos appears in both oblate 
and prolate deformations (with respect to the Kerr met- 
ric) of the MSM metric. However, he states that "it is 
impossible to present chaos regions for the oblate defor- 
mation massive bodies, because those regions found in 
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Integrator 
RK5con 
RK5var 
CCM 
CCM2 
IGEM 



e 

10~ 6 
10~ 4 
0.01 
0.01 
0.1 



result 

aborted after propagation time of r = 11254.0 because AH > 10~ e 
aborted after propagation time of r = 1087.3 because AH > 10 -6 
aborted after propagation time of r = 9559.9 because h(l, y) < 
aborted after propagation time of r = 4180.1 because AH > 10 -6 
no abortion in r = [0,50000], T ca i c = 995.3s. 



' calc 



t] until abortion 
40739.7 
149469.9 
34.2 
1353.9 



TABLE III: Results of the runs with p = 0.7 for different integration schemes. 




-0.4 



20 40 60 80 

P 

FIG. 6: The Poincare section z = for m = 2.904, a = 1.549, 
q = 0, fj, = and b = 0.8, where E = 0.971 and L z = 9.3. 
The embedded panel is a detail of the Poincare section. 

[26] are actually inside a neutron star". Even though, 
Han discussed the oblate case m = 2.904, a = 1.549, 
q = 0, /i = 0, b = 0.8, E = 0.971, L z = 9.3, but he didn't 
show the respective Poincare section. This Poincare sec- 
tion is shown in Fig. 6 and it is obvious that the main 
island of stability, which lies far from the central object, 
is embedded in a strong chaotic layer. Thus, contrary 
to Han's statement "for realistic astrophysical objects 
that are oblate" , there are chaotic orbits of "test particles 



around single oblate deformation neutron stars described 
by" the MSM model. 



6. CONCLUSION 

Applying several well established standard integration 
schemes to the system of differential equations describing 
geodesic motion in an example of a non-Kcrr spacetime 
background (i.e. MSM [24]), we showed that these inte- 
grators cannot guarantee satisfactory long-term behav- 
ior for orbits in these non-integrable systems. Therefore, 
we introduce a new integration scheme appropriate for 
evolving orbits of such systems. 

The new integration scheme effectively conserves the 
integrals of motion that such Hamiltonian systems pos- 
sess, and it is well-behaved in the case of long term evo- 
lution of strong chaotic orbits. Thus, the new integration 
scheme is well-suited for studying geodesic orbits in the 
case of spacetime backgrounds which are non-integrable 
perturbations of the Kerr spacetime. 

Moreover, we showed that chaos can appear in oblate 
deformations of the MSM metric modeling the exterior 
of a neutron star. 
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